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Summary. A driven granular material, e.g. a vibrated box full of sand, is a station- 
ary system which may be very far from equilibrium. The standard equilibrium sta- 
tistical mechanics is therefore inadequate to describe fluctuations in such a system. 
Here we present numerical and analytical results concerning energy and injected 
power fluctuations. In the first part we explain how the study of the probability 
density function (pdf ) of the fluctuations of total energy is related to the character- 
ization of velocity correlations. Two different regimes are addressed: the gas driven 
at the boundaries and the homogeneously driven gas. In a granular gas, due to non- 
Gaussianity of the velocity pdf or lack of homogeneity in hydrodynamics profiles, 
even in the absence of velocity correlations, the fluctuations of total energy are non- 
trivial and may lead to erroneous conclusions about the role of correlations. In the 
second part of the chapter we take into consideration the fluctuations of injected 
power in driven granular gas models. Recently, real and numerical experiments have 
been interpreted as evidence that the fluctuations of power injection seem to satisfy 
the Gallavotti-Cohen Fluctuation Relation. We will discuss an alternative interpre- 
tation of such results which invalidates the Gallavotti-Cohen symmetry. Moreover, 
starting from the Liouville equation and using techniques from large deviation the- 
ory, the general validity of a Fluctuation Relation for power injection in driven gran- 
ular gases is questioned. Finally a functional is defined using the Lebowitz-Spohn 
approach for Markov processes applied to the linear inelastic Boltzmann equation 
relevant to describe the motion of a tracer particle. Such a functional results to be 
different from injected power and to satisfy a Fluctuation Relation. 
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1 Introduction 

In equilibrium thermodynamics one characterizes the stable phases of a sys- 
tem using a limited set of macroscopic state variables, therefore bypassing 
much of the microscopic details of the systems under study. It is only very re- 
cently that the same strategy has been applied to systems in Non-Equilibrium 
Steady-States (NESS) . And the need for such an approach is all the more preg- 
nant for the study of NESS that no general formalism parallel to the standard 
equilibrium Gibbs-Boltzmann ensemble theory exists. This field started with 
experiments carried out on turbulent flows or convection cells, and much more 
recently on granular systems. Global observables, namely spatially integrated 
over the whole system, and their distribution, may indeed coarse-grain the ir- 
relevant microscopic details specific to the system at hand, while allowing for 
comparisons between different systems. They are expected to be more robust 
and more exportable tools for analysis than local probes, like e.g. structure 
factors. This has led to the observation of intriguing similarities between tur- 
bulent flows and granular systems [BHP98,BdSMRM05,Ber05]. However, one 
must take into account the key ingredient making a NESS way different from 
its equilibrium counterpart: steady flows of energy, matter, or else, run across 
the system. The existence of currents characterizes a NESS, and makes it dif- 
ferent from an equilibrium state in that detailed balance (time reversibility) 
no longer holds. Given that the time direction plays a central role, one is led to 
the idea that time integration may also be useful in smoothening out various 
details of the microscopic dynamics. This has motivated several authors to 
consider the distribution of time integrated and spatially averaged quantities 
characterizing the NESS as such, like that of the injected power in a turbulent 
flow or in a granular gas. 

We briefly turn to a reminder of phenomenological thermodynamics of 
nonequilibrium systems, as presented in [dGM69]. There, for systems only 
slightly away from equilibrium, the concept of entropy can be extended in a 
consistent fashion, and its time evolution goes according to 



The intrinsic entropy production rate a- m is positive definite, and cancels 
under the condition that the system reaches equilibrium. The other piece in 
the rhs of (1), which features an entropy current 3s, conveys the existence of 
external sources, often located at the system's boundaries, driving the system 
out of equilibrium. The entropy current is not but a linear combination of the 
various currents flowing through the system, with the conjugate affinities (like 
a temperature or a chemical potential gradient) as the proportionality factors. 
The entropy current - when it can univoquely be defined - therefore stands 
as a relevant measure of how far the system is from equilibrium. For that 
reason, various studies, starting from the pioneering work of Evans, Cohen 




(1) 
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and Morriss [ECM93], in a study of a thermostatted fluid under shear, have 
been focused on the appropriately generalized expression of the latter entropy 
current. In Evans, Cohen and Morriss' case it is simply proportional to the 
power provided by the thermostat to compensate for viscous dissipation. They 
went on to determine the distribution function of Qs(t), the time integrated 
entropy flow (or equivalcntly the energy provided by the thermostat), denoted 
by P(Qs,t). In doing so they empirically noticed a remarkable property of 
the pdf of Qs, namely 

1 P(Q s ,t) 

^I 1os pFq^)=^ (2) 

where qs = Qs/t, which is a time-intensive quantity, is the time average of Js 
over [0, t\. This symmetry property of the pdf of Qs was soon to be formalized 
into a theorem for thermostatted systems by Gallavotti and Cohen [GC95], 
and has since triggered a flurry of studies. The mathematical object defined 
by n(qs) = hriit-,00 log p ^ s t,f ) is seen to be extending the concept of intensive 
free energy to a nonequilibrium setting, and will occupy much of our numeri- 
cal and analytical efforts. 

In the realm of nonequilibrium systems, granular gases play a central role 
as systems exhibiting a strongly irreversible microscopic dynamics due to in- 
elastic collisions, and for these no viable definition of entropy, let alone entropy 
flow, is available. This has led various authors [AFMP01,FM04] to conjecture 
that, by analogy to thermostatted systems, the power injected into the system 
to maintain it in a steady-state, could satisfy a symmetry property like the 
one uncovered by [ECM93], and its ensuing consequences in terms of gener- 
alized fluctuation-dissipation theorems. Fortunately, a well-controlled kinetic 
theory-based statistical mechanics exists for dilute gaseous systems, and we 
shall build upon it to investigate the questions raised above. 

The outline of the present review is as follows. We begin in Sec. 2 with a 
brief introduction to granular gases and the basics of their statistical mechan- 
ics. In Sec. 3 we analyze the distribution of the total kinetic energy of the gas, 
as a first choice for a global observable. In Sec. 4 and 5 we address numeri- 
cally, analytically and also experimentally, the issue of interpreting the power 
injected into a gas in terms of entropy flow, the negative outcome of which 
leads us to Sec. 6. There we construct a one particle observable exhibiting the 
properties expected from an entropy flow, and quite notably its distribution 
function displays the symmetry property (2). 



2 A brief introduction to granular gases 



A granular gas is an assembly of macroscopic particles kept in a gaseous steady 
state by a constant excitation [BTE05] (the typical example can be illustrated 
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thinking of many beads in a strongly vibrated box). The simplest way to 
characterize such systems is to consider N identical smooth hard spheres, 
losing a part of their kinetic energy after each collision. The total momentum 
is conserved in collisions, and only the normal component of the velocity is 
affected. Thus, the collision law for a couple of particles (1,2) reads: 

v* =vi - |(1 + a) ( Via -a)* ,^ 
v* 2 =v 2 + |(1 + a)(vi2 • <t)<t , 

where <r is a unitary vector along the center of the colliding particles at 
contact. Here a is a constant, called the coefficient of normal restitution (0 < 
a < 1, and when a — 1 collisions arc purely elastic). Without some energy 
injection mechanism the total energy of the gas will decrease in time, until 
all the particles are at rest (cooling state). However, when some energy input 
is provided, the system can reach a noncquilibrium stationary state. Energy 
injection may be supplied in several ways, which can be divided in two main 
categories: injection from the boundaries and homogeneous driving. In the 
former category energy is supplied by a boundary condition, the system hence 
develops spatial gradients and it is not homogeneous. The latter category 
refers to systems where energy injection is achieved by a homogeneous and 
isotropic force acting on each particle. 

2.1 Boundary driven gases 

In this section we will give a short introduction to the methods used to describe 
the behavior of a granular gas in which energy is injected by a boundary 
condition (typically a vibrating wall). This kind of system has been widely 
studied in the literature [GZBN97, MB97, ML98, Kum98, BRMM00, BT02], 
and one of its main characteristics is that the density and the temperature 
are not homogeneous over the system: there is a heat flux, which does not 
verify Fourier law. This feature is well described by kinetic theory and in good 
agreement with the hydrodynamic approximation, which allows an analytical 
calculation of the density and temperature profiles. In the dilute limit, such a 
system is well described by the Boltzmann equation: 

a,/(r,vi,t)+vi- V/(r,vi,t) = J[f\f] . (4) 

Here J[f\f] is the collision integral, which takes into account the inelasticity 
of the particles: 

J[/|/] = cr^|dv 2 |'d^(via^) ( /(Vr '^{ (Vr,t) -/(vi ; t)/(v 2 ,^ , 

(5) 



where the notation V12 denotes the relative velocity between particles 1 and 
2, the two stars superscript (i.e. v**) denote the precollisional velocity of a 
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particle having velocity v, and the primed integral is a short-hand notation 
meaning that the integration is performed on all angles satisfying v 12 • <x > 0. 
The hydrodynamic fields arc defined as the velocity moments: 

n(r,t)= J dv/(r,v,t), (6) 

n(r,t)u(v,t) = J dvv/(r,v,t), (7) 

±n{v, t)T{T,t) = J dv|(v-u) 2 /(r,v,t), (8) 

and the hydrodynamic balance equations for those quantities are derived tak- 
ing the velocity moments in the equation (4). Their expression is: 

d t n + V • (nu) = 0, (9) 

(<9 t + u- V)ui + (mn) _1 VjPy = 0, (10) 



(d t + u • V + C) T + — {PijVjUi + V • q) = 0, (11) 



2 

3n 

where the pressure tensor Pij, heat flux q, and the cooling rate £ are defined 
by: 



(r,t)= f dvm(vi-Ui)(vj - Uj)f(r,v,t) , (12) 
q(r,i) = |dv^(v-u) 2 (v-u)/(r,v,t) , (13) 

C(r ' t} = Jr^HM)T(r!t) / dvi / dV2 ^ KW'MW)- 

(14) 

Explicit analytical expressions for the above quantities have been obtained 
in the limit of small spatial gradients by Brey et al. [BDKS98, BC01]. More- 
over for systems in the steady state without a macroscopic velocity flow the 
hydrodynamic equations simplify, and therefore the temperature and density 
profiles can be explicitly computed. 



2.2 Randomly driven gases 



We consider here a granular gas kept in a stationary state by an external ho- 
mogeneous thermostat, the so called "Stochastic thermostat", which couples 
each particle to a white noise. Energy injection is hence achieved by means of 
random forces acting independently on each particle, and drives the gas into a 
non-equilibrium steady state. The equation of motion governing the dynamics 
of each particle is therefore: 
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m^=F- u + Ff, (15) 

where F™ 11 is the force due to collisions and F' h is a Gaussian white noise 
(i.e. (F^(t)F$(t')) = 2rS ij 5^ 5 5(t - t'), where the subscripts i and j are 
used to refer to the particles, while 7 and S denote the Euclidean compo- 
nents of the random force). This model is one of the most studied in granular 
gas theory and reproduces many qualitative features of real driven inelastic 
gases [WM96, P098, PLM+98, vNETP99, HBBOO, MSS01, PTVNE02, vNE98]. 
After a few collisions per particle the system attains a non-equilibrium sta- 
tionary state. This state seems homogeneous. From the equations of motion 
it is possible to derive the homogeneous Boltzmann equation governing the 
evolution of the one-particle velocity distribution function [vNE98] : 

ft/(vi,*) = J[/,/] + /Vi Vl /(vi) , (16) 

where the Laplace operator A v = (d/d v ) 2 is a diffusion term in velocity space 
characterizing the effect of the random force, while J[/, /] is the collision inte- 
gral, which takes into account the inelasticity of the collisions (cf. eq. (5)). The 
granular temperature of the system is defined as usual as the mean kinetic 
energy per degree of freedom, T g = (v 2 )/d. The stationary solution of equa- 
tion (16) has extensively been investigated in the last years. Even if an exact 
solution is still missing, a general method is to look for solutions in the form of 
a Gaussian distribution multiplied by a series of Sonine polynomials [CC60] : 

/ st (v)=e-4 ^l + f>P 5 pQ^)) ■ (17) 

The expression of the first three Sonine polynomials is: 

S (x) = 1 

S^x) = -x + id (18) 

S 2 (x) = \x 2 - \(d + 2)x + \d{d + 2) . 
Z Z o 

Moreover the coefficients a p are found to be proportional to the averaged 
polynomial of order p: 

a P = A p( S p{^)) > ( 19 ) 

where A p is a constant and the angular brackets denote average with weight 
fst- From this observation one directly obtains that the first coefficient a\ van- 
ishes by definition of the temperature. A first approximation for the velocity 
pdf is therefore to truncate the expansion up to the second order (p = 2). An 
approximated expression for the coefficient a 2 has been found as a function 
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of the restitution coefficient a and the dimension d [vNE98, CDPT03, MSOO]. 
Its expression is: 

a ((y] = i6(I-a)(I-2 Q 2 ) 

21 ' 73 + 56d-24ad-I05a + 30(l-a)a 2 ' 1 ' 

It must be noted that the second Soninc approximation is only valid for not 
too large velocities, since the tails of the pdf have been shown [vNE98] to be 
overpopulated with respect to the Gaussian distribution. It is known [vNE98] 
that in high energies log/(i>) ~ — (v/v c ) 3 ^ 2 with a threshold velocity v c that 
diverges when the dimension d goes to infinity. This means that at high di- 
mensions the distribution is almost a Gaussian, since both the tails and the 
a,2 contributions tend to vanish. All the above results have been confirmed by 
numerical simulations, in particular through Molecular Dynamics (MD) and 
Direct simulation Monte Carlo (DSMC) [Bir94] methods. Those two numer- 
ical methods, although very different, show a surprisingly good agreement. 
This points out to correctness of the molecular chaos assumption and thus 
to the relevance of the DSMC method, which is particularly well adapted to 
simulate the dynamics of a homogeneous dilute gas. 



3 Total energy fluctuations in vibrated and driven 
granular gases 

3.1 The inhomogeneous boundary driven gas 

In this section we will study the energy fluctuations of a granular gas in 
the case where the energy is injected into the system by a vibrating wall. 
Recently Aumaitrc et al. [AFFM04] investigated, by means of Molecular Dy- 
namic (MD), the fluctuations of the total energy of the system. In particular 
they looked at the behavior of the first two moments of the energy pdf when 
the system size is changed, at constant averaged density. Because of the inho- 
mogeneities, the mean kinetic energy is no more proportional to the number 
of particles, and thus it is not an extensive quantity, and analogously the 
mean kinetic temperature is no more intensive. This has led to the definition 
of an effective (intensive) temperature and an effective number of particles, 
which makes the energy extensive. In the following we will show how a rough 
calculation (neglecting correlations and small non-Gaussianity) using the hy- 
drodynamic prediction for the temperature profile [VPB + 06b], can explain 
the phenomenology observed in [AFFM04]. Within this description it is pos- 
sible to get an expression of the effective temperature and number of particles 
as a function of the system parameters (i.e. number of particles, restitution 
coefficient, and temperature of the vibrating wall). 
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Energy Probability Distribution Function 

In this part we will compute the energy pdf for a granular gas between two 
(infinite) parallel walls. The distance between the two walls is denoted by H, 
oriented along the x axis. Here we assume that one of the walls (in x = 0) 
has small and random vibrations, acting as a thermostat that fixes to To the 
temperature at x = 0. Our boundary conditions therefore are: 



T(£ = 0) = T , 



dT 



dT 

dx 







(21) 



--H 



where the rescaled length I will be defined below. For the particular case of 
a steady state without macroscopic velocity flow, is it possible to solve those 
balance equations and get the temperature profile [BRMM00] : 



T{£) = To 



cosh I y/ a{a) (£ r , 



cosh I yj a(a)£ r] 



(22) 



where a(a) is a function of the restitution coefficient (its complete expression 
is given in rcf. [BdSMRM04]). The variable £ is proportional to the integrated 
density of the system on the x axis. Its definition is given by the following 
relation involving the local mean- free-path \{x): 



d£ 



da; 



,A(.t) 



-i -l 



d-l 



r[(d + 1)/2] 



n(x) 



(23) 



In the following we will suppose the velocity distribution to be a Maxwcllian (a 
small non-Gaussian behavior exists, but it is not relevant for this calculation) 
with a local temperature (variance) given by (22): 



e 2T («) 



{2-KT{£)) d / 2 ' 

The distribution for the energy of one particle (e = v 2 /2) is hence: 

p(e,T(£)) = f ih>§ (e) , 

where f a ,v(x) is the gamma distribution [Fel71]: 



(24) 



(25) 



r(v) 



(26) 



Our interest goes to the macroscopic fluctuations integrated over all the sys- 
tem. Thus, the macroscopic variable of interest is the granular temperature 
T g , defined here as the average of the local temperature over the x profile: 
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Tg = if L n(r)T(r) dr = h C m de ■ (27) 

with 

^ =JV - W N -=wi> (28) 

where Vd-i is the area of the surface of dimension d — 1 orthogonal to the 
x-direction, i.e. H x V^-i = V. When d = 2 one has Vd_i = L where L is 
the width of the system. N x is the number of particles per unit of section 
perpendicular to the x axis. To get an expression of the energy pdf over the 
whole system, it is useful to divide the box in £ m /A£ boxes of equal height (in 
the £ scale) A£. It is helpful to use the length scale £ because the number of 
particles Ng in each box of size L x Al is a constant. Moreover, in each box i 
we will suppose the temperature a constant Tj = T(i A€), defined expanding 
the granular temperature in a Riemann sum: 



t m /At 

- E T * M (29) 



T„ = lim 

At 



The calculation of the pdf of the box energy i.e. a sum of the energies of 
the Ni particles in a box i, is hence straightforward when the velocities of the 
particles are supposed to be uncorrelatcd: 

Qi(y) = prob(e l =y)=fj_ dNe (y) , (30) 
The characteristic function of qi (y) is 

Qi(k) = 1 —^w T ( 31 ) 

Thus, the characteristic function for the kinetic energy of the whole system 
E = ^2 £i can be obtained as the product of the characteristic functions q, (fc): 

tm/At i m /Al 

p(k)= n ^(k)= n ■ ( 32 ) 

j=q j=o (1 - ikTj) 2 

Since the number of particle per box Ng is a known fraction of the total 
number of particles (Ng = NA£/£ m ), one can rewrite the expression (32) 
as a Riemann sum. In the limit A£ — ► this yields the total kinetic energy 
characteristic function: 

P(k) = exp (-^- n [™ H (1 - i^W) A ■ (33) 

Note that this result is valid for any temperature profile T(£) and hence it 
can be applied also to other situations with different boundary conditions or 
different hydrodynamic equations. 
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Comparison with simulations 

Aumaitre et al. [AFMP01, AFFM04] showed by Molecular dynamic simula- 
tions, that the pdf of the total energy is well fitted by a x 2 l aw n{E) = 
f 1 n* (E) with a number of degrees of freedom Nf different from dN, and a 

temperature Te different from the granular temperature T g . The two param- 
eters Nf and Te are functions of the first two cumulants of the pdf: 

N '^wt- T ' = W' (34) 

The notation (X) c denotes the cumulant of the variable X. Here we want to 
compare result (33) with these numerical results. Since we are not able to 
analytically calculate the inverse Fourier Transform of (33) using (22) as a 
temperature profile, we used a numerical computation to obtain it in an ap- 
proximate form. Moreover, an expression of the cumulants of the total kinetic 
energy can be obtained from the characteristic function (33): 

= w \ TP{i)de ■ (35) 

Jo 

In figure 1 the Inverse Fourier Transform of (33) is compared with the function 
n{E) previously defined. The similarity of the two functions is remarkable. 
Another important feature that can be checked with this results is the de- 
pendence of the above defined two macroscopic quantities (Nf and Te) with 
system size. It is straightforward to see, from (27) and (35), that the granular 
temperature and the total kinetic energy arc respectively an intensive and an 
extensive variable if £ m is independent from the system size. This is effectively 
the case if both the density p = N/V and the total height H are kept con- 
stant. Moreover, for large enough £ m , the integral in equation (35) becomes 
size independent: 

/ " T"(£)d£ ~ T f . (36) 
Jo 2p^a{a) 

Thus, the effective temperature Te defined above becomes a constant propor- 
tional to the temperature of the wall, while the parameter Nf still depends 
on the system size: 

1 dN T 
ya(a) *m * 

Numerical simulations show that Te effectively remains a constant for large 
systems, and under several procedures of box size increase. The behavior of 
Nf is determined by the maximum of the integrated density £ m . For a square 
cell at constant density one finds £ m oc s/N , so that Nf oc s/N, which is not 
far from N° observed in [AFFM04]. Moreover, if only the height H of the 
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cell is increased, £ m is proportional to N, and Nf becomes constant. All those 
features are in agreement with the numerical observations in [AFFM04] . The 
above results clearly show that a rough calculation, which takes into account 
only the inhomogencitics of the system, is able to quantitatively describe the 
behavior of the fluctuations of the total kinetic energy of a vibrated granular 
gas. In some cases the energy pdf can be approximated with a gamma distri- 
bution, which is the standard distribution for the energy pdf in the canonical 
equilibrium. Nevertheless there are strong deviations from the equilibrium the- 
ory of fluctuations, since the two parameters of the gamma distribution (i.e. 
the temperature and the number of degrees of freedom) are not the granular 
temperature neither the number of degrees of freedom. Another important re- 
mark is that correlations, and in particular contributions from the two points 
distribution function, do not play a primary role to explain those deviations 
from the equilibrium theory of fluctuations. In order to characterize correc- 
tions arising from the two particles velocity pdf, one should measure energy 
fluctuations at a given height x from the vibrating wall. As already noted 
in [AFFM04] this task is very hard, since the available statistic become very 
poor. Nevertheless an effective way to quantify those fluctuations is to look at 
homogeneous systems, where contributions coming from the inhomogencitics 
vanish. With this objective in mind, we will be interested in the following in 
granular gases heated by an homogeneous and isotropic driving. 



Fig. 1. Energy pdf (solid line) and a gamma distribution with same mean and 
same variance (dotted line) for a restitution coefficient a = 0.9 , N — 100 particles 
in two dimensions in a box of density p = 0.04, height H = 50, and with a wall 
temperature To = 5. 



3.2 The homogeneously driven case 

In this section we will present some numerical results concerning the energy 
fluctuations in a dilute gas driven by the stochastic thermostat presented 
in section 2.2 [VPB + 06b]. When the system reaches a stationary state, the 
dissipated energy is compensated by the energy injected by the thermostat, 
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and the temperature fluctuates around its mean value. Here we are interested 
in the fluctuations of the total energy measured by the quantity 

° E = N — JEW 2 — ■ ( } 

Note that a% = 2N/Nf. Brcy et al. have computed, by means of kinetic 
equations, an analytical expression for a\ in the homogeneous cooling state, 
which is equivalent to the so-called Gaussian deterministic thermostat. One 
of the main differences of this stochastic thermostat with a deterministic one, 
is found in the elastic limit. On the one hand, for the cooling state, when the 
restitution coefficient tends to 1, the conservation of energy imposes that the 
energy pdf is a Delta function, and the quantity &e goes to 0. On the other 
hand, with the stochastic thermostat, if the elastic limit is taken keeping the 
temperature constant, the strength of the white noise will tend to zero, but 
it will still play a role in the velocity correlation function. 

We performed DSMC simulations to measure the energy pdf of such a sys- 
tem. A plot of this pdf is shown in figure 2, and it is close to a ^-distribution 
with same mean and same variance. Nevertheless the number of degrees of 
freedom of this ^-distribution is lower than the true number of degrees of 
freedom (i.e. (N — 1) x d). This effect may arise from two separated causes: 
the non-gaussianity of the velocity pdf, and the presence of correlations be- 
tween the velocities. This feature also suggests that a calculation of the energy 
pdf with the hypothesis of uncorrelated velocities (but non-Gaussian) could 
explain at least a part of this non-trivial effect. In order to quantify these 
contributions we will consider that the velocity pdf is well described by a 
Gaussian multiplied by the second Sonine polynomial: 



f ^=i^fW{ 1 + a ^{2T)) ' (39) 

where 02 is given by expression (20). 

The calculation of the pdf of the sum of the square of N variables dis- 
tributed following (39) is straightforward. The characteristic function of the 
energy pdf is: 

~ „, 1 / d(d + 2) ( 1 2 xx 



(40) 



where N is the number of particles of the system. This yields: 



(E) = ^NT, {E 2 )-(E) 2 = ^NT 2 (l + ^a 2 ^ . (41) 

It is now possible to have an explicit expression for the energy fluctuations: 
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In figure 3 this result is compared with the result of DSMC simulations, per- 
formed for several values of the restitution coefficient a and for two different 
values for the number of particles N . The disagreement between the uncorre- 
cted calculation and the simulations is a clear sign of the correlations induced 
by the inelasticity of the system. One can note that the fluctuations increase 
when the restitution coefficient decreases. One can also see that there is a 
value of the restitution coefficient a around l/y/2, that is when the approxi- 
mate expression of a 2 vanishes, for which a\ is exactly 1 = 2/d, as for a gas 
in the canonical equilibrium (velocities are then uncorrelated) . 

We now turn to the dependence of a% on the strength of the white noise 
r. It is useful, for this purpose, to introduce a rescaled, dimensionless energy 

E = t E -W . (43) 

We have plotted in figure 4 this rescaled energy pdf for a system of N = 
100 particles with a restitution coefficient a = 0.5 and for several values of 
the strength of the white noise r. One can see how all the pdfs collapse 
into a unique distribution. The role of the noise's strength is thus only to 
set the temperature (or mean kinetic energy) scale. Besides, relative energy 
fluctuations depend only on a and N. Moreover, since a\ does not depend 
on the number of particles N (for N large enough) , the central limit theroem 
applies, and hence P(E) is a Gaussian in the thermodynamic (N — > 00) limit. 
In conclusion we have shown that randomly driven granular gases display non 
trivial fluctuations, because of the correlations induced by the inelasticity. Two 
different kinds of correlations contribute to this behavior of the fluctuations. 
First, the non-Gaussianity of the velocity pdf, which simply tells that the 
Euclidean components of the velocity of each particle are correlated one to 
each other. Second, a contribution from the two particles velocity pdf, which 
does not factorize exactly as a product of two one-particle distributions. It 
must be pointed out, however, that these correlations do not invalidate the 
Boltzmann equation. As already noted in [EC81, BdSMRM04, CPM07], the 
two points correlation function <72(vi,V2), which is defined by: 

5 2 (vi,v 2 ) = / (2) (vx,v 2 ) - /(vx)/(v 2 ) , (44) 

where is the two points distribution, is of higher order in the density ex- 
pansion (roughly speaking C(p(vi,v 2 )) ~ O (/(vi)/(v 2 )) /N). This is con- 
firmed by the numerical observations, since when the number of particles 
increases, the energy pdf tends be closer and closer to a Gaussian. Spatial 
correlations, which can be at work in homogeneously driven granular gases, 
at higher densities, and which have been negleted here (assuming spatial ho- 
mogeneity) can also play a relevant role in fluctuations [PBV07] . 
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Fig. 2. Energy pdf (dots) from DSMC simulations with a restitution coefficient 
a = 0.5 and ./V = 100 particles for a system driven with the stochastic thermostat. 
The solid line shows a gamma distribution with same mean and same variance. 
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Fig. 3. Plot of a% versus the restitution coefficient a for N = 100 (O) an d N = 1000 
(□) particles. The result of the calculation assuming uncorrelated velocities (42) is 
shown by the dashed line. 
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Fig. 4. Plot of the pdf of the rescaled energy E for a restitution coefficient a — 0.5 
and for N = 100 for several values of the strength of the noise F. 
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4 A large deviation theory for the injected power 
fluctuations in the homogeneous driven granular gas 

From now on we turn our attention to the fluctuations of another global quan- 
tity, i.e. the power injected into the system by the external source of energy. 
In particular in this section our main goal is to obtain a kinetic equation able 
to describe the behavior of the large deviations function of the time integrated 
injected power [VPB+05, VPB + 06a] in a randomly driven gas (cf section 2.2). 
The latter quantity is the total work VV provided by the thermostat over a 
time interval [0, t]: 



W(t)= f dt^Ff . v< 



(45) 



Our interest goes to the distribution of W(i), denoted by P(W, t), and to its 
associated large deviation function Ttoo(w) defined for the reduced variable 
w = W/t (W(t) being extensive in time): 

TToo H = lim TT t (w) , TT t (w) = - log P( W = wt, t) . (46) 
t— too t 

We introduce p(-T/v, W, i) the probability that the system is in state Jjv at 
time t with W(t) = W. The function we want to calculate is 



P(w,t) =1 dr NP (r N ,w,t). (47) 

We shall focus on the generating function of the phase space density 

p( r N, A, t) = J dWe~ xvv p(r N ,W, t) (48) 
and on the large deviation function of 

P(x,t)= [ dWc~ xw p(w,t) = { dr N p(r N ,\,t) (49) 



which we define as 

fj,(X) = lim -logP(A,<) . (50) 
Note that p{X) is the generating function of the cumulants of W, namely 

(51) 



lim = (_!)« *m 



A=0 



Moreover 7Too(w) can be obtained from fi(X) by means of a Legendre transform, 
i.e. TToo(w) = jU(A«) + A»w with A* such that jjf(X*) = —w. 

The observable W is non-stationary but it is Markovian, hence a general- 
ized Liouville equation for the extended phase-space density p(P/v,W,i) can 
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be written. It varies in time under the combined effect of the inelastic collisions 
(which do not alter W) and of the random kicks: 

dtp = d t p +d tP (52) 



collisions 



kicks 



Considering that the thermostat acts independently on each particle, it can 
be shown that 

dtp = V \r(A Vi + 2A/v t ■ d Vi + r{d\ + a\ 2 )1 p (53) 

kicks — * L J 

i 

This additional piece is linear in p just as the collision part is. The large 
time behavior of p is governed by the largest eigenvalue p,{X) of the evolution 
operator of p. In the large time limit, we thus expect that 

p(r N , A, t) ~ C(X)e^ t p(r N , A), (54) 

where /5(-T/v,A) is the eigenfunction associated to p, and C(A) is such that 
p(r^, X) is normalized to unity. We then introduce 

/ W (»! , . . . , v k , A, t) = J dr N - k p, (55) 

where J d/jv-ft means an integration over N — k particles, we have that 

<9 t / (1) (v, x,t) = rA v f + 2xra v ■ vf + r{x 2 v 2 - dx)f + J (56) 

with J — J dWc~ xw J the Laplace transform of the collision integral in which 
f(v, W, t) now plays the role of the velocity distribution. Quite unexpectedly 
the above equation has a straight physical interpretation: consider a many 
particle system where a noise of strength _T and a viscous friction-like force 
F = — 2ArV act independently on each particle, and where the particles inter- 
act by inelastic collisions. Consider then that the particles annihilate/branch 
(depending on the sign of A) at constant rate dXr, and branch with a rate 
proportional to X 2 v 2 T. Then, the equation governing the evolution of the one 
particle velocity distribution of such a system is exactly the equation (56), 
where A is a parameter tuning the strength of the external fields. Moreover, 
in spite of there being no a priori reason for that, p, as well as / = J d-T/v-iP, 
can be interpreted as probability density functions. 

The one and two-point functions f^'(v, W, t) and f^ 2 '(v-i,V2,W,t) that 
enter the expression of J are expected to verify, at large times, 

fW(v 1 ,\,t) = C(\)e*fW{v 1 ,\), (57) 

and 

P\v l ,v i ,X,t) = C{X)ef tt f^{v 1 ,V2,X), (58) 

where both and f^ 1 are normalized to unity. We perform the following 
molecular-chaos- like assumption: 
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f^(v 1 ,v 2 ,X)^P\v 1 ,X)f^(v 2 ,X) (59) 

which does have a definite physical interpretation in the language of the inelas- 
tic hard-spheres with fictitious dynamics (viscous friction, velocity dependent 
branching/annihilation) described in the above paragraph. Then we get that 



fj,f(v, A) =rA v f + 2APv • d v f + r{d\ + AV)/ 

1 

+ i 



tk>2d<7Vi2 • <T 

V12 -CT>0 



a' 2 fK*, A)/(«r, A) - />!, A)/> 2 , A) 

(60) 



where we have now omitted the superscript (1) denoting the one-point func- 
tion. The A = limiting case yields the usual Boltzmann equation, since in 
this case a stationary solution exists, and hence /Lt(A = 0) = 0. The boundary 
condition to the evolution equation above is thus: 

/(v,A = 0) = / st (v) (61) 

with /st(v) the stationary velocity pdf (cf. Eq.(17)). 



4.1 The cumulants 

Here we find an approximated expression of /i(A) solving a system of equations 
obtained projecting (60) on the first velocity moments. First we shall define 
a dimensionless velocity c = v/vo(X), where ^o(A) plays the role of a thermal 
velocity: 

v 2 (A) = 2T(A) = \ [ dv v 2 f(v, A) . (62) 



Then, defining the function /(c, A) = t>o(A)/(v, A), and its related moments 
of order n 

m„(A) = J dc c n f(c,X) , (63) 
one obtains the following recursion relation: 

r 

(/t + r(2n + d)X)m„ = —*n(n + d — 2)m„_ 2 + rA 2 i>nm„ + 2 — v v n7 (64) 
where 

r dcc" J[f,f] . (65) 

Recalling the definition of the cumulants (51), and the approximated solution 
for the stationary velocity pdf, it appears natural to argue that, for A ~ 0, 
the function /(c, A) should be well approximated by: 



/(c, A) = 0(c) (1 + a 1 (X)S 1 (c 2 ) + a 2 (X)S 2 (c 2 )) + 0(a 3 ) , (66) 
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where (f>(c) = 7r~ d//2 cxp(— c 2 ) is the Gaussian distribution. Even in this case, 
from the relation (19) and from the definition (62), the coefficient a\ is found 
to be 0. The method consists in taking the equation (64) for n = 0, 2 and 4 
in order to find an explicit expression of /i, «o, and 02 in the limit A — ► 0. The 
quantities 1/2 and ^4 have been calculated at the first order in 02 [vNE98] , and 
their explicit expressions are: 

(1-a 2 ) Q d f 3 1 dr f 3 1 



and 



with 



Vi = JE= {Ti + a 2 T 2 } , (68) 



T l= d+l + a 2 (69) 

r 2 = |,(10d + 39 + 10a 2 ) + ^4, (70) 
32 (1 — a) 

( \ 2/3 

where To = ( w~^^f A ) ' s the the granular temperature obtained averag- 
ing over Gaussian velocity pdfs (i.e. the zero-th order of Sonine expansion). 
The expression of the first moments m n is: 

m = 1 (71a) 
m 2 = d/2 (71b) 

m 4 = (71c) 

m 6 = (l + So^t + d) < 4 + d > (71d) 
8 

With the help of the above defined temperature scale To, we introduce some 
dimcnsionless variables: 

A* = Mttt , A = AT , 

(72) 



2T ' T 



Note that this scaling naturally defines the scales for the other quantities of 
interest, namely: 

T n w ~ W 

= * = dT< W = W)- (73) 

The expression of the moment equation (64) becomes, for the above defined 
dimcnsionless quantities: 
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fid + (2n + d)XJ m n = ^ m n -2 + 2v m n +2 - vov n . (74) 

First we solve the above equation for n = 0, getting the following result: 

/}(A) = -A + A 2 S 2 (A). (75) 

Recalling that when A — > one has Vq = 2T g + O(X), it is important to note 
that if we restrict our analysis to the Gaussian approximation for P(W,t), 
that is if we truncate /j,(X) to order A 2 , eq. (75) will read: 

^ = A(AT 3 - 1) . (76) 



Then we see that indeed 



fi(A) = m (-r - , (77) 

which means that 7roo(u>) = maxA{^(A) + Xw} verifies 

w 

7Too(w) - TToo^—w) = — . (78) 

1 a 

However, the nontrivial functions m n (X) will break the property (77), as we 
shall explicitly show later. In order to characterize more precisely the depen- 
dence of jj, upon A for small values of A , it is useful to expand Vq and a2 in 
powers of A: 



5g(A) =«5 +\v£ +X z v£ +0{X' i ) (79a) 
02(A) = 4 0) + a4 1} + A 2 4 2) + 0(X 3 ) (79b) 

In this way we can find w 2< ' ( a 2^ solving equation (74) for 



= 2: 



~2 (0 > 
^0 



2 d) _ 4 4 0) 4 1} 



3 



(80) 
(81) 



, r >. 2 _ af (i. + !) + f _f . (82) 

Then we substitute Wq(A) in the third equation and expand it in powers of A 
to find the expression of 0,% (a) . Note that one has also to expand in powers 
of <i2 and keep only the linear terms in order to be coherent with the v p 
calculations. We find the following expressions, which arc plotted in Fig. 5: 
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,(°) 



4 (1 - a) (1 -2 a 2 ) 



19 + 14 d - 3 a (9 + 2d) + 6(1 - a) a 2 

4(1 -a) 2 (-l + 2a 2 ) (31 + 2a 2 + 16d) 
(19 + 14 d - 3 a (9 + 2d) + 6(1 - a) a 2 ) 2 



(2) _ A{a) 
B(a) 



a 2 = 



(83) 
(84) 
(85) 



with 



A(a) = 16 (-1 + a) 2 (-1 + 2 a 2 ) {906 + a [-984 + a (85 + 3 a (-19 + 6 (-1 + a) a))] + 985 d+ 

+ a [-951 + a (-25 + 3 a (7 + 6 (-1 + a) a))] d + (269 + 3 a (-75 + 2 a (-7 + 3 a))) d 2 } , 

(86) 

and 

B(a) = 3 (-19- 14d + 3a (9 + 2 (-1 + a) a + 2d)f (87) 

The v 2 } ' expression, as well as the expression, gives the usual results 
established for granular gases [vNE98, MS00]. At this point the computation 
of the cumulants becomes straightforward. From relation (51) it follows: 



= (-i) n jvdny _1 n! £ 2< "~ 2) . 



lim 

t^oo t 

Moreover, since the corrections are numerically small, the zcro-th order 
(Gaussian) approximation already gives a good estimate for the cumulants. 
Namely, the first cumulants arc , in this approximation: 

(W) r = tNdr , (W 2 ) = 2tNdrT , 

(89) 

(W 3 ) c = 8tNdrT$ , ( w4 ) c = 48iiVdrT 3 . 

All the above expansions in powers of A, at the second order in Sonine co- 
efficients (e.g. 02) can be carried out just expanding vq and 02 in (79) to 
higher powers of A. Moreover, expanding in higher order in Sonine coefficient 
(e.g. 03) remains in principle still possible, but it will involve a higher number 
of equations in the hierarchy (74) (e.g. n = 6), and therefore will need the 
expression of higher order collisional moments (e.g. v§). 



4.2 The solvable infinite dimension limit 

Strong arguments [VPB + 06a] can be given showing that in high dimensions 
f(v,X) is not far from a Gaussian. We are therefore led to consider, in the 
limit d — ► 00, /(v, A) to be a Gaussian with a A-dependent second moment. 
In this situation the dimensionless function / will read: 
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4-th order 
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Fig. 5. 

d = 2. 



ai 1 ' and a, versus q for 



Fig. 6. The solid line shows fi in 
the limit d — > oo. The dashed line is 
/t at fourth order in A from (75) for 
d = 2 and a = 0.5. Finally the dot- 
ted line shows the same quantity cal- 
culated with a truncation at second or- 
der in A, which would satisfy the G-C 
relation. 



/(c) 



n d/2 



(90) 



with c = v/vo(A). In this context one can solve equation (74) in order to 
get an explicit expression for ^t(A). Solving the system defined by (74) for 
n = and n = 2 gives a unique solution for /2(A) which verifies the physical 
requirement /i(0) = 0: 



with: 



So(A) = 



1 + 4 A 3 6i(A) 
2 A 4 



32 
A2 



-A + y5 2 (A). 



(l + 4A 3 )' 



and 



6i(A) = A /A" 8 + ^ - b 2 (\) + 6 3 (A) , 63(A) 



+ 62(A) - 63(A) + 



is (§r 



(91) 



64(A) 
46i(A) 

(92) 



A 5 



63(A) = 



2(1) 



9 + V3 V27 + 256A3 



A 4 



6 4 (A) = ™ 
V ; A 3 



A 3 ( 9 + V3 V27 + 256A 3 



192 ( 1 + 4 A 3 



1+4A 3 



A ( > 



(93) 

This expression of the velocity scale reduces to the kinetic temperature for 
A = 0, and decreases monotonically as A -1 / 2 when A — ► 00. This means that 
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in the limit A — > oo / approaches a Dirac distribution as cxp(— Xv 2 /2). This 
feature supports the intuition that the small W events (which are related to 
the large values of A) are provided by the small velocities. The behavior of 
ft is shown in Fig. 6. The large deviations function ft(X) becomes complex 

for A < — ;p73 , because of the terms containing y27 + 256A 3 . Moreover for 
large A the behavior of this function is ft(X) ~ —A3. In the vicinity of the 
singularity (i.e. A = Ao = — ^ e behavior of the large deviation function 
is: 

£W = - 3 2/3 2 1/6 \/A - Ao + 0{\ - A ) . (94) 

From the behavior for large A it is possible to recover the left tail of the large 
deviation function -Koq. In general, if /i(A) ~ — \P for A — > oo, this leads to 
//(A*) = — /3A» _1 = —w. This last relation tells us that for /3 < 1 we are 
recovering the limit w — > + , with a behavior of the large deviation function 
given by ttoo(w) — /i(A«) + \*w ~ u" 3 - 1 . Moreover, from the behavior of /i 
near Ao, an analogous calculation provides the right tail of the large deviation 
function: ■n 00 (w) ~ Aow, when w —> oo. Finally, in our particular case, the 
tails are given by 

7Too(w — * + ) ~ -W~ 1/3 , TTooiw — > Oo) ~ -tZ) , (95) 

Note that there is no w < tail to ir^. The graph of the whole function 
ti'oo('U') is depicted in Fig. 7. 

o 
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5 Fluctuations of injected power at finite times: two 
examples 

5.1 The homogeneous driven gas of inelastic hard disks 

In this section the results of numerical simulations of two models (inelastic 
hard spheres and inelastic Maxwell model) are presented with particular at- 
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tcntion to the verification of the Fluctuation Relation for the injected power. 
The main requirement to pose the question about the validity of the Fluc- 
tuation Relation is a clean observation of a negative tail in the pdf of the 
injected power. This dramatically limits the time t of integration of W(t). In 
numerical simulations, as well as in real experiments, at time larger than a 
few mean free times the negative tail disappears. On the other hand, at times 
of the order of 1-3 mean free times, the Fluctuation Relation appears to be 
correctly verified for the inelastic Hard Spheres model and slightly violated 
for the inelastic Maxwell model. The measure of the cumulants, anyway, gives 
a neat indication of the fact that the time of convergence of the large devia- 
tion function is at least 10 times as large and that the true asymptotic is well 
reproduced by the theory exposed in this chapter. This theory shows strong 
arguments against the validity of a symmetry relation of the Gallavotti-Cohen 
type for the large deviations of injected power. 

The stationary state of a driven granular gas, modeled by equation (16), 
under the assumption of Molecular Chaos may be studied with a Direct Sim- 
ulation Monte Carlo technique [Bir94,MS00]. As a first check of reliability of 
the algorithm, we have measured the granular temperature T g and the first 
non-zero Sonine coefficient a-i = ((v 4 ) / (v 2 ) 2 — 3)/3. The measured granular 
temperature is always in perfect agreement with the estimate. The measured 
a-2 coefficient is a highly fluctuating quantity and its average is in very good 
agreement with the theoretical estimate. 

In figure 8 the probability density functions p(w, i) = tP(wt, t) (for t equal 
to 1 mean free time) for three different choices of parameters N, r (at fixed 
restitution coefficient a) is shown. The values of the first two cumulants of the 
distribution and their theoretical values are compared in table 1, with very 
good agreement. In the same table we present also the measure of the third 
and fourth cumulants. 



N 


r 


{W(t))/t 


<w(t) 2 )cA 


NTd 


2NrdT g 


Wt) 3 )c/( 


<W(i) 4 >cA 


100 


0.5 


100 


20835 


100 


21052 


6.02779 x 10 b 


1.54181 x 10 s 


100 


12.5 


2500 


13019125 


2500 


13157900 


9.47684 x 10 a 


6.12963 x 10 la 


200 


0.5 


199.9 


42009 


200 


42120 


1.21911 x 10 b 


3.09634 x 10 s 



Table 1. Rescaled cumulants of the distribution of injected work P(W, i), measured 
with t equal to 1 mean free time for different choices of the parameters. 



The comparison with a Gaussian with same mean value and same variance 
shows that the pdf P(W, t) is not exactly a Gaussian. In particular there are 
deviations from the Gaussian form in the right (positive) tail. This is well seen 
in figure 9. It must be noted that the important deviations in the right tail 
arise at values of W{t) larger than the minimum W(t) available in the left 
tail, i.e. they have no influence in the following plot of figure 10 regarding the 
Gallavotti-Cohen symmetry. 
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Fig. 8. Probability density func- 
tion of the injected power, p(w,t) = 
tP(W(t) = wt,t) with t equal to 1 
mean free time. In all three cases the 
value of the restitution coefficient is 
a = 0.9. Other parameters are a) TV = 
100, r = 0.5; b) N=100, r = 12.5; 
c) N = 200, r = 0.5. The dashed 
line represents a Gaussian with same 
first two cumulants. These distribu- 
tions have been obtained with ~ 1.5 x 
10 9 independent values of W(t). 




Fig. 9. Ratio of P(W, t) and a Gaus- 
sian with the same first two moments, 
for the same parameters as in figure 8: 

a) corresponds to N — 100, r = 0.5, 

b) to TV = 100, r = 12.5 and c) to N = 
200, r = 0.5. The range between the 
vertical dotted lines is the useful one 
for the check of the Gallavotti-Cohen 
relation. It can be noted that the 
strongest deviations from the Gaussian 
behavior appear outside of this range. 



In figure 10 the Gallavotti-Cohen relation iTt(w) — n t (—w) — fleffw is 
questioned for the same choice of the parameters. The relation, at this level of 
resolution and for this value of the time t (1 mean free time), is well satisfied. 
Moreover table 2 shows that the value of (3 e f / is well approximated by /3 = 
1/Tg, as expected if the truncation of /x(A) at the second order were valid, 
see eq. (78). In figure 11 the same relation is checked for different values 
of t, slightly larger (i.e. up to t equal to 3 mean free times). No relevant 
deviations are observed as t is increased. Moreover this figure is important to 
understand the dramatic consequences that a larger t has on the "visibility" 
of the Gallavotti-Cohen symmetry: as t is increased, events with negative 
integrated power injection become rarer and rarer. This eventually leads to 
the vanishing of the left branch of P(W,t). 



N 


r 


Peff 


l/T g 


100 


0.5 


0.0100 


0.00955 


100 


12.5 


0.000402 


0.000382 


200 


0.5 


0.00995 


0.00952 



Table 2. Factor of proportionality in the "Gallavotti-Cohen" relation compared 
with j3. 
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w 



Fig. 10. Plot of nt(w) — 7rt(— w) vsw, 
for the same choice of the parameters 
as in figure 8. The values of the slope 
peff of the fitting dashed lines are in 
table 2 




Fig. 11. Plot of nt(w) — 7Tt(— w) ver- 
sus w, for the system with N = 100 
and r = 0.5 for different values of t. 
We recall that in this case (w) = 100. 
The dashed line has slope (3 = 1/T g . 
In the inset the corresponding p(w, t) 
are shown. 



The main conclusion is that no appreciable departure from the A 2 trunca- 
tion is observed at this level of resolution. Much larger statistics are required 
to probe the very high energy tails of p(w, t). Further numerical insights make 
evident that the small times used to check the GC Relation (t smaller or equal 
than 3 mean free times) are far from the time where the asymptotic large de- 
viation scaling starts working. In figure 12 we show indeed the numerical 
measure of the third cumulant of W(t) rescaled by the first cumulant, varying 
the integration time t. The time of saturation is of the order of ~ 50 mean 
free times. The saturation value is in very good agreement with the value 
predicted by our theory, eq. (89). Note that this value is not at all trivial, 
since the third cumulant for a Gaussian distribution is zero. At that time the 
measurable 7r t (w) is shown in figure 13, rescaled by (w). The accessible range 
of values from a numerical simulation is dramatically poor and we think it is 
already remarkable to have obtained a good measure of the third cumulant 
with such a resolution. 

The reason for a verification at small times of the GC formula is the 
following: near w = the pdf of w is almost a Gaussian. In the Gaussian case 
wc immediately get Kt(w) — 7it(— w) = f3 e ffW with f} e ff — 2(W(t)) / (W{t) 2 ) c . 
The first two cumulants at small times are easily obtained considering an 
uncorrelated sequence of energy injection, obtaining (W(t))/t = NTd and 
(W(t) 2 ) c /t = F-' 1 • v ?; ) 2 ) c = 2NTdT g . Then the value p eff = 1/T g is 
unavoidable. In this case the GCFR observed is nothing else than the Green- 
Kubo (or Einstein) relation, which is known to be valid for driven granular 
gases: (W(t) 2 )/t = 2T g (W(t))/t [PBL02, PBV07]. Small deviations from a 
Gaussian appear, in first approximation, as small deviations from the slope 
1/T g , but the straight line behavior is robust since the first non-linear term 
of ir t (w) - ir t (-w) is not w 2 but w 3 [AFMP01]. 
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Fig. 12. Third cumulant W = 
W/((W)T 9 2 ) for several times of inte- 
gration. The time is in units of the 
mean collision time. Note that the time 
when a stationary value of the rescaled 
cumulant is reached is much larger 
than the characteristic time of the sys- 
tem (the collision time). 
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Fig. 13. Numerical measurement of 
7Tt for a time of 50 collisions per par- 
ticle (when a stationary value for the 
rescaled cumulant is reached). 



Numerical simulations of the Inelastic Maxwell Model have been performed 
with a Direct Simulation Monte Carlo analogous to the one used in the Hard 
Spheres model. The Maxwell gas is a kinetic model due to Maxwell, who ob- 
served that a pair potential proportional to r - 2 ( d - 1 ) j r being the distance 
between two interacting particles, gives rise to a great simplification of the 
collision integral [Max67]. In fact this kind of interaction makes the collision 
frequency velocity independent. It must be noted that when the inelasticity 
of the particles is considered, this model looses its straight physical inter- 
pretation, but it nevertheless keeps its own interest. The collision integral is 
analytically simpler than the hard particles model and preserves the essential 
physical ingredients in order to have qualitatively the same phenomenology. 
In the recent development of granular gases this kinetic model has been ex- 
tensively investigated [BMP02,BNK02,BNK03,EB02,BCG00]. Thanks to the 
simplifications present in this model, we are able to improve the number of 
collected data by more than a factor of ten. The distributions of the injected 
power p(w, t) are shown in figure 14 for some choices of the restitution coeffi- 
cient a. The driving amplitude r has been changed in order to keep constant 
the stationary granular temperature T g . In figure 15 we have displayed the 
deviations from the Gaussian of P(W,t). The non-Gaussianity of P(W, t) is 
highly pronounced, but again it is striking only in the positive branch of the 
pdf. We have tried, with success, a fit with a fourth order polynomial, which 
is consistent with the usual truncation of the Soninc expansion to the second 
Soninc polynomial. 

Finally, in figure 16, we have attempted a check of the Gallavotti-Cohen 
fluctuation relation. The relation seems to be systematically violated. This 
appears in two points: 1) the right-left ratio of the large deviation function 



Fluctuations in granular gases 27 



is not a straight line; 2) the best fitting line has a slope which is larger than 
l/Tg. The "curvature" (and the deviation from the 1/T g line) increases with 
decreasing values of a, indicating that the inelasticity is the cause of the devi- 
ation from the Gallavotti-Cohen relation. It should be noted that to achieve 
this result we have collected more than 4 x 10 10 independent values of W(t), so 
that the statistics of the negative large deviations could be clearly displayed. 
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Fig. 14. p(w,t) = tP(wt,t) for dif- 
ferent values of a (at fixed constant 
temperature T g ) in the Driven Inelas- 
tic Maxwell Model measured at a time 
t equal to 1 mean free time. The 
dashed lines are Gaussian distributions 
with the same mean and same vari- 
ance. These distributions have been 
obtained with ~4x 10 10 independent 
values of W(t). 
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Fig. 15. p{w, t) (at t equal to 1 mean 
free time) divided by a Gaussian with 
same average and same variance for 
different values of a (at fixed constant 
temperature T g ) in the Driven Inelas- 
tic Maxwell Model. The light dashed 
lines represent a fit with a polynomial 
of fourth order. 



5.2 The boundary driven gas of inelastic hard disks 

In a recent experiment on vibrated granular gases [FM04] it has been ar- 
gued that the statistics of the power injected on a subsystem by the rest of 
the gas fulfills the Fluctuation Relation (FR) by Gallavotti and Cohen. The 
experiment was performed by putting in a two-dimensional vertical box TV 
disks of glass and submitting the container to a strong vertical vibration. We 
have reproduced the experiment by means of a Molecular Dynamics (MD) 
simulations of inelastic hard disks, observing perfect agreement with the ex- 
perimental results and obtaining a deeper insight into the system. The main 
difference of this model with respect to the previous "homogeneously driven" 
model is that the external energy source is located at the two horizontal (top 
and bottom) walls. This boundary driving mechanism leads to the develop- 
ment of spatial inhomogeneities and the appearance of internal currents. 

The event driven MD simulations have been performed for a system of TV 
inelastic hard disks with restitution coefficient a, diameter a and mass m = 1. 
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Fig. 16. (Color online). Finite time check of Gallavotti-Cohen relation for the 
injected power (with t equal to 1 mean free time), i.e. 7r t (iy) — -Kt(-w) vs. w, in 
a numerical simulation of the Driven Inelastic Maxwell Model with N = 50, and 
different values of a (the driving amplitude r has been rescaled in order to fix the 
granular temperature T g ). The dashed curve is a straight line with slope /3 — 1/T g . 
The dotted curve is a straight line obtained fitting the a = 0.1 data points until 
w = 45, useful as a guide for the eye. The thin (red) solid curve is a fit with a cubic 
(0.28w + 5.6 ■ 10~ i w 2 - 1.1 ■ 10" 5 w 3 ). 
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Fig. 17. Left: Snapshot of the system considered for MD simulations, with the 
inner region marked by the solid rectangle. Right: Corresponding vertical profiles 
of density (&(y), dashed line) and temperature (T(y), solid line). The dotted lines 
mark the bottom and top boundaries of the inner region. Here iV = 270 and a = 0.9. 
The mean free path is ~ 5.7d. 
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The vertical 2D box of width L x = 48a and height L y = 32a is shaken by a 
sinusoidal vibration with frequency / (period Tb ox = 1//) and amplitude 2.6a. 
Collisions with the elastic walls inject energy and allow the system to reach 
a stationary state. We have checked that possible inelastic collisions with the 
walls hardly affect the results. Gravity - set to g = —1.7af 2 in order to be 
consistent with the experiment- has a negligible influence on the measured 
quantities. We have varied the restitution coefficient from 0.8 up to 0.99 (glass 
beads yield on average a w 0.9) and the total area coverage from 0.f38 (i.e. 
N = 270) up to 0.32 (N = 620). In figure 17-left a snapshot of the system is 
shown. During the simulations the main physical observables are statistically 
stationary. The local area coverage field ^(x, y) and the granular temperature 
field T(x,y) (defined in 2D as the local average kinetic energy per particle) are 
almost uniform in the horizontal direction, apart from small layers near the 
side walls. In figure 17-right the profiles <P(y) — (1/L X ) J dx<l > (x, y) and T(y) = 
(1/L X ) J dxT(x, y) are shown to be symmetric with respect to the bottom and 
the top of the box. Following the experimental procedure, we have focused our 
attention on a "window" in the center of the box, fixed in the laboratory frame, 
of width 2L x /5 and height L y /3, marked in figure 17-left. Apart from the 
negligible change of potential energy due to gravity, the total kinetic energy 
of the particles inside the window, changes during a time r because of two 
contributions: AK T = Q T — I T where Q T is the kinetic energy transported 
by particles through the boundary of the window (summed when going-in 
and subtracted when going-out) and I T is the kinetic energy dissipated in 
inelastic collisions during time r. For several values of r we have measured, 
as in the experiments, Q T which is related to the kinetic contribution to the 
heat flux (we checked that inclusion of the collisional contribution, even if 
non small [HMZ04], docs not change the picture). With N = 270 and a = 0.9 
the characteristic times are the mean free time r co ; ~ 0.47t6 o;c , the diffusion 
time across the window Tdiff = 0.82ti, ox and the mean time between two 
subsequent crossings of particles (from outside to inside) T cross w 0.039n ox . 

We define the injected power as q T = Q t /t and two relevant probabil- 
ity density functions (pdfs): /q(Q t ) an d fq(q T )- figure 18a) shows f q {q T ) for 
different values of r. A direct comparison with Fig. 3 of Ref. [FM04] sug- 
gests a fair agreement between simulations of inelastic hard disks and the 
experiment. The pdfs are strongly non-Gaussian and asymmetric, becoming 
narrower as r is increased. At small r a strong peak in q T = is visible. 
More interestingly, f q {q T ) at small values of r has two different exponen- 
tial tails, i.e. f q (q T ) ~ exp(=F/3-tT<jv) when q T — > ±oo with /?_ > /3+. The 
peak and the exponential tails at small r are observed also in the experiment 
(see Fig. 3 of [FM04]) and in similar simulations [AFMP01]. In figure 18b) 
we display log[f q (q T )/f q (—q T )]/T vs. q T , which is equivalent to the graph of 
tTt{<1t) — Tr T (—q T ) vs. q T where Tr T (q T ) = log[/Q(rg T )]/r. From figure 18 it ap- 
pears that at large values of r, 7r r (g T ) — 7r T (— q T ) is linear with a r-independent 
slope (3 e ff 7^ 1. We have measured /3 e // with various choices of the resti- 
tution coefficient a and of the covered area fraction finding similar results. 
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Fig. 18. a) pdfs of injected power f q {q T ) from MD simulations for different values 
of r = (1,2,4,8,16,32) x T mi „ with T min = 0.015tw Here N = 270 and a = 0.9. 
The distributions are shifted vertically for clarity. The dashed lines put in evidence 
the exponential tails of the pdf at r = Tmin- b) plot of (1/r) log[/g(g-r)//<z(— <fr)] 
vs. q T from MD simulations (same parameters as above) at large values of r. The 
solid curve is a linear fit (with slope f3 e ff) of the data at r = 128r m i„. The dashed 
line has a slope /3 gra n = 1/T gran . In the inset the same graph is shown for small 
values of r = (1, 2, 4, 8) x Train (from bottom to top). 

Ref. [FM04] reports f3 e ffT gran ~ 0.25 where T gran is the mean granular tem- 
perature in the observation window. Similar values are measured in our MD 
simulations. At area fraction 13.8% and a = 0.9 we have P e ffT gra n ~ 0.23. 
At fixed a and increasing area fraction, f3 e ffT gran slightly increases, as in the 
experiment. As a — ► 1 the slope Peff decreases. At a = 1 (without gravity and 
external driving) the distribution of Q T is symmetrical and (3 e ff = 0, indicat- 
ing that 1 / p e f / is not a physically relevant temperature concept. Interestingly, 
it appears that (3 e ff is a non hydrodynamic quantity: different systems may 
show the same density and temperature profiles, with very different values of 



We now adopt a coarse-grained description of the experiment which is able 
to entirely capture the observed phenomenology. The measured flow of energy 
is given by 



where n_ {n+) is the number of particles leaving (entering) the window dur- 
ing the time r, and vf_ (vf + ) are the squared moduli of their velocities. In 
order to analyze the statistics of Q T we take n_ and n + as Poisson-distributed 
random variables with average lot, where to corresponds to the inverse of the 
crossing time T cross . In doing so we neglect correlations among particles en- 
tering or leaving successively the central region. A key point, supported by 
direct observation in the numerical experiment, lies in the assumption that the 
velocities Vi+ and come from populations with different temperatures 7+ 
and T_ respectively. Indeed, compared with the population entering the cen- 



Peff- 




(96) 
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tral region, those particles that leave it have suffered on average more inelastic 
collisions, so that T_ < T + , Finally we assume Gaussian velocity pdfs. Within 
such a framework, the distribution /q(<5 t ) of Q T can be studied analytically. 
Here it is enough to recall that i v i > m D dimensions, if each compo- 

nent of Vi is independently Gaussian-distributed with zero mean and vari- 
ance T, is a stochastic variable x with a distribution Xti,t(x) = fi/T,Dn/2(x), 
where f a ,v{x) is the Gamma distribution, and whose generating function reads 
Xn,r(z) = (1 — Tz)~ Dn ^ 2 [Fcl71]. It is then straightforward to obtain the gen- 
erating function of Q T in the form Jq(z) = exp[r^(z)] with 

H(z) = w (-2 + (1 - T+z)- D ' 2 + (1 + T_z)- D / 2 ) . (97) 

We observe that Jq(z) has two poles in z = ±1/T± and two branch cuts on 
the real axis for z > 1/7+ and z < — l/TL. From fi(z) we immediately obtain 
the cumulants of Jq(Q t ) through the formula (Q n ) c = t ^"M(0)- 

For r — ► 00 the large deviation theory states that Jq(Q t ) ~ exp(r7r 00 ((5 T /T)) 
and Tr 00 (q) can be obtained from /i(z) through a Legendre transform, i.e. 
""00(9) = max(fj,(z) — qz). The study of the singularities of /j,(z) reveals the 
behavior of the large deviation function iv 00 (q) for q — > ±00. It can be seen 
that 

q q 
ttooO?) — (g->oo), TTooiq)^— (q^-00). (98) 

We emphasize however that it is almost impossible to appreciate these tails 
in simulations and in experiments, since the statistics for large values of q and 
t is very poor. 

A Gallavotti-Cohen-type relation [GC95,LS99,Kur98], e.g. ^^(q)— 7r 00 (— q) = 
(3q for any q and an arbitrary value of (3 would imply fj,(z) = /x(/3 — z). 
One can see that such a (3 does not exist, i.e. the fluctuations of Q T do not 
satisfy a Gallavotti-Cohen-like relation. The observed linearity of the graph 
log[fq(q T )/f q (—q T )]/T = n(q T ) — Tr(—q T ) vs. q T can be explained by the same 
observation pointed out in the previous subsection: at large values of r it is 
extremely difficult, in simulations as well as in experiments, to reach large 
values of q, while for small q, Tr(q) — tt(— q) w 2w'(0)q + o(q 3 ), i.e. a straight 
line with a slope j3 e ff = 2V (0) is likely to be observed. It has been already 
shown [Far02] that in dissipative systems deviations from the FR can be hid- 
den by insufficient statistics at high values of q. The knowledge of fi(z) is 
useful to predict this slope. At large r, 7r'(0) w J7'(0) = — z*(0) where z*(q) 
is the value of z for which n(z) — qz is extremal. This gives 

When T + = XI (i.e. if a = 1) /3 e // = 0. We emphasize that /3 e // does not 
depend upon w. We have compared with success these predictions with the 
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numerical and experimental results, measuring the temperatures T + and T_ 
in the simulation. 

What happens for small values of r? We note that /q(z) has the form 
cxp (r/i ( z )) for any value of r and not only for large r. Therefore at small r one 
can expand the exponential, obtaining /q(z) ~ 1+wr (—2 + (1 — T + z)~ D / 2 + 
(1 + T_z)~ D / 2 ) . This immediately leads to an analytical expression for fq{Q T ) = 
const x 8{Q T ) + Xi,t + (Qr) + Xi,t_ (— Qt), which fairly accounts for the strong 
peak which is observed in the experiment and in the simulations, and predicts 
exponential tails for /q(Q t ): Xi,t{x) oc x d I 2 ~ x exp(— x/T) so that (3 + = 1/7+ 
and /3_ = 1/T_. This suggests an experimental test of this theoretical ap- 
proach: the measure at small values of r of the slopes of the exponential tails 
of /q(Qt) should coincide with a direct measure of T+ and T_. However, we 
point out that the values of /3+ and /?_ obtained by fitting the tails in the 
hard disks simulation, using values as small as r = 0.00015ti, ox yield esti- 
mates of 7+ and T_ which are larger (by a factor ~ 1.6) than those found 
by a direct measure. This disagreement brings the limits of such a simple 
two-temperature picture to the fore. In the simulation and in the original ex- 
periment the measured injected energy is indeed the sum of several different 
contributions, namely Q T w Q% x + Q% v + Q^ x + Q yy where Q l J is the kinetic 
energy transported by the i component of the velocity by particles crossing the 
boundary through a wall perpendicular to direction j. Two main differences 
with the simplified interpretation given above arise: a) there are two couples 
of temperatures, i.e. T|,Tf as well as T\,T V _ [BC98,ML98,BK03]; b) the di- 
agonal contributions Q J T J are sums of squares of velocities whose distribution 
is not a Gaussian but is ~ v exp(— v 2 /T), since the probability of crossing is 
biased by the velocity itself. The calculation of /q(Q t ) is still feasible, with 
qualitatively similar results. 

6 The dynamics of a tracer particle as a non-equilibrium 
Markov process 

In the search for a quantity that is, more rigorously, related to the "en- 
tropy production" in a granular gas, we consider in this section the pro- 
jection of the dynamics of the gas onto that of a tracer particle, which is 
easier since it is equivalent to a jump Markov process. We are interested 
in the dynamics of a tracer granular particle in a homogeneous and dilute 
gas of grains which is driven by an unspecified energy source. The require- 
ments are that the gas is dilute, spatially homogeneous and time transla- 
tional invariant. The gas is characterized by its velocity probability den- 
sity function P(v) which, for the sake of simplicity, will be considered of 

the form P(v) = (27r j )d/2 exp + a 2 S${v 2 /2T)), (where S$ is the 

d-dimcnsional second Sonine polynomial already defined in 2.2). The gas is 
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therefore parametrized by its temperature T and its second Sonine coefficient 
a 2 which measures its non-Gaussianity. 

The linear Boltzmann equation for the tracer, in generic dimension d, 
reads: 



dP*(v,t) 1 



j J &V X / dV2 j d ^l( V l - V 2) • U>|P*(vi)P(v 2 )x 

x h (v - Vl + ^^[(vi - v a ) • - J(v - vi) 



(If 



where P*(v) is the velocity pdf of the test particle and the primed integral 
again indicates that the integration is performed on all angles that satisfy 
(vi — V2) • u> > 0. The mean free path I appears in front of the collision 
integrals. In the following (when not stated differently) we will put 1=1, 
which can be always obtained by a rcscaling of time. 

We rewrite the above equation (100) as a Master equation for a Markov 
jump process [PVTVW06]: 



dP.(v,t) 
df 



= y"dv 1 P >t (v 1 )A"(v 1 ,v)-y"dv 1 P,(v)A"(v,v 1 ). (101) 



The transition rate K(y, v') of jumping from v to v' is given by the following 
formula: 



K{w,w') = (j^) \\M 2 - d J dv 2T P[v 2 (v,v',v 2T )], (102) 

where Av = v' — v denotes the change of velocity of the test particle after a 
collision. The vectorial function v 2 is defined as 

v 2 (v, v', v 2t ) = v 2(7 (v,v')<t(v,v') + v 2t , (103) 

where <x(v, v') is the unitary vector parallel to Av, while v 2t is entirely con- 
tained in the (d— l)-dimcnsional space perpendicular to Av (i.e. v 2r -ZW = 0). 
This implies that the integral in expression (102) is (d — l)-dimensional. Fi- 
nally, to fully determine the transition rate (102), the expression of u 2( j is 
needed: 

u 2CT (v,v') = —— \Av\ + v & . (104) 
1 + a 

6.1 Detailed balance 

Here, we obtain a simple expression for the ratio between if(v,v') and 
K{y' , v). When exchanging v with v' the unitary vector <r changes sign. Fur- 
thermore one has that w 2(T (v,v') 7^ v 2o -(v',v). From all these considerations 
and from equation (102) one obtains immediately: 
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AT(v,v') |dv 2r P[v 2 (v,v')] _ PKM] 

tf(v',v) Jdv 2T P[v 2 (v',v)] -PK(v',v)] ' 1 U ° j 

We note that this ratio depends only on the choice of the pdf of the gas, P, 
and not on the other parameters (such as a). However in realistic situations 
(experiments or Molecular Dynamics simulations) P is not a free parameter 
but is determined by the choice of the setup (e.g. external driving, material 
details, geometry of the container, etc.). 

Introducing the short-hand notation v 2o - = «2cr(v,v'), v' 2a = «2cr(v',v) 
and i£ = • <x, we also note that 

«f = 4a + K + v' a f - 2n 2a {v a + v' a ) , (106) 
from which it follows that 

A 2 = (v 2a f «) 2 = -A- 2 l —^A = , (107) 

1 + a 1 + a 

where A = v 2 — {v' a ) 2 = \v\ 2 — \v'\ 2 , i.e. the kinetic energy lost by the test- 
particle during one collision. When a = 1 then A 2 — —A (energy conserva- 
tion). From the above considerations it follows that 

• in the Gaussian case, it is found 

, if(v.v') A I- a A 3- a A ,„ . 

log— H — r = — + 2 = 108) 

6 A"(v',v) 2T l + a2T \ + a2T v ' 

• in the First Sonine Correction case, it is found 



log 



K(v,v') _ 3- a A 
K(v',v) ~ 1 + aXf 



log 



1 + a 2 S$ 



1 



,TT^( v ' < t- v ")+ v ") 2 
2T 



IT 



(109) 



In the case where P(v) is a Gaussian with temperature T, it is immediate 
to observe that 

P*(v)/v(v,v') = P*(v')#(v',v) (110) 

if P» is equal to a Gaussian with temperature T" = < T. This means 

that there is a Gaussian stationary solution of equation (101) (in the Gaussian- 
bulk case), which satisfies detailed balance. The fact that such a Gaussian 
with a different temperature T" is an exact stationary solution was known 
from [MP99]. It thus turns out that detailed balance is satisfied, even out of 
thermal equilibrium. Of course this is an artifact of such a model: it is highly 
unrealistic that a granular gas yields a Gaussian velocity pdf. As soon as 
the gas velocity pdf P(v) ceases to be Gaussian, detailed balance is violated, 
i.e. the stationary process performed by the tracer particle is no more in 
equilibrium within the thermostatting gas. We will see in section 6.2 how to 
characterize this departure from equilibrium. 
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6.2 Action functionals 

From the previous section we have learnt that the dynamics of the velocity 
of a tracer particle immersed in a granular gas is equivalent to a Markov 
process with well defined transition rates. This means that the velocity of 
the tracer particle stays in a state v for a random time t > distributed 
with the law r(v)e~ r ( v ' t dt and then makes a transition to a new value v' 
with a probability r(v) -1 A(v, v'), with r(v) = J dv'K(v,v'). At this point 
it is interesting to ask about some characterization of the non-equilibrium 
dynamics, i.e. of the violation of detailed balance, which we know to happen 
whenever the surrounding granular gas has a non-Gaussian distribution of 
velocity. 

To this extent, we define two different action functionals, following [LS99]: 

p(vi -> v 2 -> ... -> v„ (t) ) 
P(v n (t) -> v n ( t )_i -> ... -> vi) 

where i is the index of collision suffered by the tagged particle, is the velocity 
of the particle before the i-th. collision, v' { is its post-collisional velocity, n(t) 
is the total number of collisions in the trajectory from time up to time 
t, and K is the transition rate of the jump due to the collision. Finally, we 
have used the notation P(v± —> V2 v„) to identify the probability 

of observing the trajectory vi — * v 2 — > ... — * v„. The quantities W(t) and 
W(t) arc different for each different trajectory (i.e. sequence of jumps) of 
the tagged particle. Note that the first term log p P ^7 1 ' > j in the definition of 

Wit), cq. (111c), is non-extensive in time. The two above functionals have the 
following properties: 

• W(i) = if there is exact symmetry, i.e. if K(vi — > v%+i) = A'(v.; + i — > v^) 
(e.g. in the microcanonical ensemble); W(t) = if there is detailed balance 
(e.g. any equilibrium ensemble); 

• we expect that, for large enough t, for almost all the trajectories lim^oo W(s) / s 
linis-xx, W(s)/s = (W(t)/t) = (W(i)/t); here (since the system under in- 
vestigation is ergodic and stationary) the meaning of () is intuitively an 
average over many independent segments of a single very long trajectory; 

• for large enough t: 1) at equilibrium (W(t)) = (W(t)) — 0; 2) out of 
equilibrium (i.e. if detailed balance is not satisfied) those two averages are 
positive; we use those equivalent averages, at large t, to characterize the 
distance from equilibrium of the stationary system; 
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• if S(t) = — J dvP*(v, t) log P*(i>, t) is the entropy associated to the pdf of 
the velocity of the tagged particle P*(v, t) at time t (e.g. —H where H is 
the Boltzmann-H function), then 

±S(t)=R(t)-A(t) (112) 

where R(t) is always non-negative, A(t) is linear with respect to P* and, 
finally, (W(t)) = f Q At' Alt'). This leads to consider W(t) equivalent to the 
contribution of a single trajectory to the total entropy flux. In a stationary 
state A(t) = R(t) and therefore the flux is equivalent to the production; 
this property has been recognized in [LS99] . 

• FRw (Lebowitz-Spohn-Gallavotti-Cohen fluctuation relation): 7r(u;) — 
7r(— w) = w where tt(w) = limt^oo j log fw(tw) and f w (x) is the probabil- 
ity density function of finding W(t) ~ x at time t; at equilibrium the FRw 
has no content; note that in principle Tr'(w,t) = j log fw(tw) 7^ ir(w) 
at any finite time; a generic derivation of this property has been ob- 
tained in [LS99] , while a rigorous proof with more restrictive hypothesis is 
in [Mae99] ; the discussion for the case of a Langevin equation is in [Kur98] . 

• FRyy- (Evans-Scarles fluctuation relation): W(w,t) —lf(—w,t) = w where 
W(w,t) = i log f^y(tw) and f^-{x) is the probability density function of 
finding W(t) = x at time t; at equilibrium the FR^ has no content; 
this relation is derived in [LS99]; the analogy between this relation and 
the Evans-Searles fluctuation relation [ES94, ES02] has been put forward 
in [PVTVW06] 

A detailed numerical study [PVTVW06] of the fluctuations of W{t) and 
W(t) in this model has shown on the one hand that, out of equilibrium (i.e. 
when the surrounding gas is non-Gaussian), the FRyy is always satisfied. 
On the other hand the FRw is always violated, even if it was expected on 
the basis of the arguments given in [LS99]. The difference between the two 
functionals defined in (111) is a term which is non-extensive in time, but 
which has fluctuations whose distribution has exponential tails and therefore, 
in principle, can contribute to the large deviation function of W(t) [PRV06, 
Vis06]. Such a failure of a large time Fluctuation Relation, which is much 
more pronounced in the near-to-equilibrium cases, is similar to that observed 
in other systems [ESR05, Far02, vZC03, BGGZ05] 



7 Conclusions 

The study of the fluctuations of global physical quantities in a granular 
gas is at its very beginning. In the lack of a general rigorous theory in 
the framework of 11011 equilibrium statistical mechanics, experiments and nu- 
merical simulations are the main source of results, together with few exact 
analytical calculations. In this review of recent results [VPB+05, PVB+05, 
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VPB + 06a, PVTvW06, VPB + 06b] we have indicated some routes that have 
been followed, focusing on two global quantities (total energy and energy 
injection rate) that are of interest in nowadays physics of non-equilibrium 
systems [BHP98, BdSMRM05, AFMP01 , AFFM04, Far02, Far04] . On one hand 
we have shown that total energy fluctuations have a pdf that strongly depends 
on the model considered. We have also pointed out that definitive inferences 
about the presence of correlations, starting from the observation of "anoma- 
lous" pdfs of total energy, must be drawn with caution, since the lack of 
spatial or temporal translational invariance may play a major role. On the 
other hand we have presented a method to calculate the large deviation func- 
tion of injected power in a granular gas: this method strongly suggests the 
disappearance of a negative branch in such large deviation function. This re- 
sult is a direct consequence of the time-irreversibility of inelastic collisions: 
injected power fluctuations are dominated at large times by the energy dissi- 
pated in collisions, which is always positive. Finally we have sketched a recipe 
to obtain a quantity related to time-reversal asymmetry (i.e. violation of de- 
tailed balance) whose large deviations can be both positive and negative. This 
quantity has the advantage of being measurable in experiments, but the dis- 
advantage of not having an obvious "macroscopic" counterpart. It contains in 
fact information on the non-Gaussianity of the velocity pdf of the gas. This 
consideration is in our opinion the main open issue in this study. 
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